function [config_sph, en] = change_config_sph(config_sph0, dr, dtheta, dphi)
    N = length(config_sph0);
    d = 0.88;
    config_sph = config_sph0;
    config = config_sph0;
    r_min = 1.1; r_max = 2.3;
    theta_min = 0; theta_max = pi;
    phi_min = 0; phi_max = 2*pi;
    
    for i=14:N
        pos_old = config_sph0{i};
        while 1
            %如果产生的新位置超出空间边界，采用周期边条件的方式处理
            r_new = pos_old(1) + dr * rand;
            if(r_new > r_max)
                r_new = r_new - (r_max - r_min);
            end
            theta_new = pos_old(2) + dtheta * rand;
            if(theta_new > theta_max)
                theta_new = theta_new - (theta_max - theta_min);
            end
            phi_new = pos_old(3) + dphi * rand;
            if(phi_new > phi_max)
                phi_new = phi_new - (phi_max - phi_min);
            end
            x = r_new * sin(theta_new) * cos(phi_new);
            y = r_new * sin(theta_new) * sin(phi_new);
            z = r_new * cos(theta_new); 
            
            new = [x,y,z];
            
            judge = 1;

            for j=1:(i-1)
                rj_vector = config{j} - new;
                rj_square = sum(rj_vector.^2);
                rj = sqrt(rj_square);
                if(rj < d)
                    judge = 0;
                    break
                end
            end

            if(judge == 1)
                config_sph{i} = [r_new, theta_new, phi_new];
                config{i} = new;
                break
            else

            end
        end
    end
    en = sim_en_cal(config);
    return
end

%在知道N=13时的能量的基础上，计算体系能量
function en = sim_en_cal(config)
    N = length(config);
    en = -11.08158890;
    for m=14:N
        for k=1:(m-1)
            en = en + LJ_pot(config{m},config{k});
        end
    end
    return
end

%取原子直径为0.88sigma，截断半径为3.0sigma
%pt1,pt2均为1*3的数组，记录原子中心在三维直角坐标系下的3个坐标分量
function pot = LJ_pot(pt1, pt2)
    %根据pt1,pt2两点的坐标计算两点距离
    rvector = pt1 - pt2;
    r_square = sum(rvector.^2);
    r = sqrt(r_square);
    pot = r^(-12)-r^(-6); 
    
%     if(r < 0.88)
%         pot = 0.88^(-12)-0.88^(-6);
%     elseif(r > 3)
%         pot = 0;
%     else
%         pot = r^(-12)-r^(-6);   
%     end

    return
end